Shock-Induced  Turbulence  and  Acoustic  Loading  on 
Aerospace  Structures 


Dimitris  Drikakis 


CRANFIELD  UNIVERSITY 
CRANFIELD  DEFENCE  AND  SECURITY 
COLLEGE  ROAD 

BEDFORD  MK430AL  UNITED  KINGDOM 


EOARD  Grant  #FA9550-14-l-0224 
Report  Date:  August  2015 
Final  Report  from  1  August  2014  to  31  July  2015 


Distribution  Statement  A:  Approved  for  public  release  distribution  is  unlimited. 


Air  Force  Research  Laboratory 
Air  Force  Office  of  Scientific  Research 
European  Office  of  Aerospace  Research  and  Development 
Unit  4515,  APO  AE  09421-4515 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply 
with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


Research  Project 


Shock-Induced  Turbulence  and  Acoustic  Loading  on 

Aerospace  Structures 


First  Year  Report 


High-Order  Methods  for  Supersonic  and 
Hypersonic  Shock  Wave  Turbulent  Boundary  Layer  Interaction* 


Dimitris  Drikakis 
University  of  Strathclyde, 
Glasgow,  UK 

dimitr  i  s .  dr  ikaki  s  @  strath  .ac.uk 


An  investigation  of  the  accuracy  of  high-order  methods  for  hypersonic  shock  wave 
turbulent  boundary  layer  interaction  (SWTBLI)  is  presented.  The  numerical  methods 
considered  here  comprise  of  the  Monotone-Upstream  Central  Scheme  for  Conservation 
Laws  (MUSCL)  and  Weighted  Essentially  Non-Oscillatory  (WENO)  schemes,  2"^  to  9^*^  order 
accurate  in  conjunction  with  structured  and  mixed  element  unstructured  grids.  Both 
Implicit  Large  Eddy  Simulation  (ILES)  and  Reynolds  Averaged  Navier-Stokes  (RANS) 
computations  have  been  performed.  The  effects  of  discretization  on  the  turbulence  transport 
equation,  including  the  approximation  method  for  the  viscous  gradient,  are  also  investigated. 
The  accuracy  of  the  schemes  in  high  Reynolds  number  RANS  modeling  is  assessed  against 
experimental  data  of  a  shock  impingement  on  a  flat  plate  at  Mach  number  5  and  unit 
Reynolds  number  37xl0Vm.  ILES  has  been  performed  for  the  compression  ramp  case  at 
moderate  Reynolds  numbers  of  38.7x10^  based  on  the  boundary  layer  thickness,  and 
compared  to  Direct  Numerical  Simulations  (DNS). 


I.  Introduction 

SHOCK-WAVE  turbulent  boundary  layer  interaction  (SWTBLI)  is  of  particular  interest  to  structural  engineers 
for  the  design  and  manufacture  of  aerospace  structures.  Pulsating  flows  featuring  unsteadiness  attributed  to 
SWTBLI  can  lead  to  fatigue  and  structural  damages\  Advancing  our  understanding  of  SWTBLI  and  associated 
loading  is  important  for  developing  effective  control  strategies  that  will  mitigate  these  loads. 

Numerical  simulations  of  SWBLI  flows  are  constrained  by  accuracy  and  computational  cost.  The  accuracy  is 
hampered  by  excessive  numerical  dissipation  and  dispersion  errors.  First  and  second  order  numerical  schemes  are 
highly  dissipative,  thus  leading  to  incorrect  predictions  of  turbulent  SWTBLI  induced  separation.  High-order  (HO) 
schemes  in  conjunction  with  large  eddy  simulations  lead  to  significantly  better  results,  however,  pressure 
fluctuations  can  be  under-,  or  over-predicted,  due  to  numerical  dispersion  errors.  In  the  last  few  years,  however, 
significant  progress  has  been  made  with  regards  to  the  application  of  high-resolution  (HR)  and  HO  methods  to 
compressible  flows  featuring  acoustic  excitation,  turbulent  SWTBLI,  and  low-Mach  number  effects. 


^  The  work  presented  in  this  report  was  carried  out  during  the  first  year  of  the  project,  while  the  PI  was  at  Cranfield 
University,  Cranfield,  MK43  OAL,  UK 
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The  simulations  and  results  presented  in  this  report  concern  Cases  1  and  2  of  the  original  program  of 
work,  i.e.  a  shock  of  Mach  number  5  (Re=37xl0^)  impinging  on  a  flat  plate  agitating  the  turbulent  boundary  layer 
and  forming  a  separation  zone;  and  implicit  large  eddy  simulation  of  supersonic  flow  over  a  compression  ramp.  In 
the  past,  several  efforts  have  been  made  to  obtain  reliable  measurements  for  supersonic  boundary  layers^.  The 
reported  experimental  results  encompass  an  alarming  degree  of  scatter  for  (nominally)  compatible  measurements 
performed  in  different  facilities.  Schiilein^  performed  several  experimental  activities  of  SWTBLI  flow  at  Mach  5. 
The  state-of-the-art  techniques  were  employed  to  measure  wall  pressure  loadings,  skin  friction,  velocity  profiles  and 
heat  loads.  This  case  has  been  extensively  employed  for  validation  and  investigation  of  computational  methods.  The 
computed  pressure  loads  are  usually  in  good  agreement  with  the  experiment,  however,  considerable  uncertainty  is 
found  in  the  skin  friction  predictions  downstream  of  the  boundary  layer  separation.  In  this  paper,  2"^  and  3*^^  order 
MUSCL,  and  3^^  and  5^^  order  WENO  numerical  schemes,  are  employed  to  compute  the  supersonic  turbulent  flow  at 
two  shock-generator  angles:  10  and  14  degrees.  The  RANS  and  ILES  results  are  compared  with  experimental  data^ 
and  DNS^^ 


II.  Governing  Equations 


The  governing  equations  are  the  3D  compressible  Navier-Stokes  equations,  which  can  be  written  in  the  following 
Cartesian  coordinates  form  after  neglecting  external  body  forces: 

W  is  the  vector  of  the  conserved  variables;  and  Fc  and  are  the  inviscid  and  viscous  fluxes,  respectively: 
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where  p  is  the  density,  u,  v  and  w  are  the  velocity  components,  p  is  the  pressure  and  E  is  the  total  energy  per  unit 
mass.  The  contravariant  velocity  is  given  by: 
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where  Xy  are  the  shear  stresses  and  k  is  the  heat  conductivity  constant.  An  ideal  gas  is  assumed  for  the 
thermodynamic  closure  of  the  equations  and  the  Sutherland’s  Eaw  is  used  in  the  calculation  of  the  dynamic 
viscosity. 
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III.  Numerical  Methods 


A.  Compressible  unstructured  grid  framework 

The  equations  are  discretized  by  a  k-exact  Finite  Volume  (FV)  method  on  mixed-element  unstructured  meshes. 
The  MUSCL  and  WENO  schemes  are  employed  for  reconstructing  the  element  averaged  solutions  using  high-order 
stencils  and  polynomial  functions.  A  detailed  description  of  the  methods  and  their  implementation  can  be  found 
in^’^.  A  brief  description  of  the  methods  is  presented  below. 

The  numerical  approach  adopted  in  the  present  study^^’^^  is  suitable  for  unstructured  meshes  with  various  types 
of  element  shapes  in  2D  and  3D.  It  has  been  previously  used  successfully  for  laminar,  transitional  and  turbulent 
flows^’^^.  A  Gaussian  numerical  quadrature  of  appropriate  order  for  the  order  of  the  polynomial  used  is  implemented 
for  the  approximation  of  the  integral  expressions  of  the  fluxes.  The  calculation  of  the  numerical  convective  and 
viscous  fluxes  requires  the  knowledge  of  the  pointwise  values  of  the  conserved  vector  as  well  as  of  the  velocity  and 
temperature  gradients  at  each  Gaussian  integration  point.  These  pointwise  values  are  approximated  through  an 
interpolation  (reconstruction)  procedure  of  a  desired  order  of  accuracy  utilizing  the  cell  averages.  The  latter  requires 
a  recursive  stencil  construction  process  where  the  direct  side  neighbor  elements  are  added  until  a  target  number  M  of 
stencil  elements  has  been  reached.  For  MUSCL-type  schemes  only  one  central  stencil  is  used,  while  for  WENO 
schemes,  in  addition  to  the  central  stencil,  several  additional  directional  stencils  are  also  employed. 

The  reconstruction  is  carried  out  in  a  transformed  system  of  coordinates  in  order  to  minimize  scaling  effects 
that  appear  in  stencils  consisting  of  elements  of  different  sizes,  as  well  as  to  improve  the  condition  number  of  the 
system  of  equations^^’^^.  For  computing  the  degrees  of  freedom,  a  minimum  of  K  cells  are  needed  in  the  stencil 
in  addition  to  the  target  cell.  Using  the  minimum  possible  number  of  cells  in  the  stencil  (M  =  K)  has  been  found  to 
produce  ill-conditioned  systems^^’^^'^^,  hence  the  choice  to  use  M  =  2K  improves  the  robustness  of  the  method. 
This  is  especially  worthwhile  since  no  substantial  performance  penalty  is  incurred  as  a  result  of  this 
improvement^^’^^’^\  The  resulting  least-squares  system  is  solved  by  a  QR  decomposition  and  the  reconstruction 
polynomial  is  computed. 

In  the  present  study  two  different  schemes  are  employed  for  discretizing  the  convective  fluxes  of  the  equations:  a 
MUSCE  scheme  using  a  TVD-type  slope  limiter^^  and  different  WENO  schemes^^’^^  based  on  the  characteristic 
variables.  For  the  viscous  part  a  linear  reconstruction  polynomial  of  the  same  order  for  the  velocity  and  temperature 
field  is  constructed  using  the  same  central  stencil  as  for  the  conserved  vector.  The  discontinuous  states  of  the 
convective  fluxes  are  approximated  by  the  HEEC  Riemann  solver^^,  and  the  central  averaging  approach  is  used  for 
the  discontinuous  viscous  flux.  The  solution  is  advanced  in  time  by  an  explicit  TVD  Runge-Kutta  B’^^-order 
method^"^.  It  is  worth  mentioning  that  unstructured  grids  for  complex  geometries  can  benefit  when  combined  with 
variational  optimization  techniques'^  as  well  as  with  very  high  order  methods^"^.  For  RANS  computations,  the 
Spalart-Allmaras  (SA)  one-equation  turbulence  model  has  been  employed. 

B.  Compressible  structured-grid  framework 

The  structured-grid  code  solves  the  full  Navier- Stokes  equations  using  a  finite  volume  Godunov-type  method. 
The  inter-cell  numerical  fluxes  are  computed  by  solving  the  Riemann  problem  using  the  reconstructed  values  of  the 
conservative  variables  at  the  cell  interfaces.  The  reconstruction  stencil  is  a  one  dimensional  swept  unidirectional 
stencil  (ID  split).  The  Riemann  problem  is  solved  using  the  Harten,  Eax,  van  Eeer,  and  (the  missing)  “Contacf’ 
(HEEC)  approximate  Riemann  solver The  reconstructed  values  utilized  in  the  HEEC  Riemann  solver  are  obtained 
primarily  by  two  different  limiter  approaches,  the  Monotone  Upstream-centered  Schemes  for  Conservation  Eaws^"^ 
(MUSCE)  and  the  Weighted-Essentially-Non-Oscillatory^^  (WENO)  reconstruction  methods.  The  MUSCE  is 
employed  with  accuracy  up  to  5^^  order^^’^^  (henceforth  labeled  M5),  whereas  the  WENO  schemes  uses  up  to  9^^- 
order  of  accuracy^^’^^  (henceforth  labeled  W9).  The  9^^-order  WENO  is  implemented  in  conjunction  with  a  relative 
smoothness  limiter^^  as  well  as  a  mapping  technique^\ 

All  the  reconstruction  techniques  can  be  further  augmented  with  a  low-Mach  limiting  scheme^,  which  involves 
an  additional  stage  in  the  reconstruction  process  for  the  velocity  vector.  This  low-Mach  number  correction  (labeled 
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as  LM)  ensures  uniform  dissipation  of  kinetic  energy  in  the  limit  of  zero  Mach  number,  thus  extending  the  validity 
of  Godunov-type  methods  to  at  least  Mach-lO""^  via  a  progressive  central  differencing  of  the  velocity  components.  It 
has  been  recently  shown^^’^^  that  use  of  the  low-Mach  number  correction  can  lead  to  a  minor  reduction  in  accuracy 
when  used  along  with  the  9*-order  WENO.  It  was  found  that  when  the  numerical  dissipation  is  sufficiently  small, 
dispersive  errors  originating  from  the  WENO  reconstruction  can  become  dominant.  Nonetheless,  the  9^^-order 
WENO  was  also  shown  to  possess  a  remarkable  inherent  low-Mach  number  capability  as  it  was  able  to  successfully 
resolve  flow  features  at  relatively  low  Mach  number  regimes  (Mach~0.1)  at  moderate  grid  resolutions  without  the 
use  of  any  low-Mach  number  correction  method.  In  the  light  of  the  above,  in  the  present  study  the  low-Mach 
number  correction  method^  was  implemented  only  in  conjunction  with  the  5^^-order  MUSCE  scheme  (henceforth 
labeled  M5EM). 

The  viscous  part  of  the  equations  is  solved  using  a  second  order  central  difference  scheme.  Finally,  the  solution 
is  advanced  in  time  using  a  three-stage  total  variation  diminishing  (TVD)  Runge-Kutta  (RK)  method^"^. 


IV.  Results 


A.  Impinging  shock 

The  case  is  based  on  the  experiment  of  Schiilein^  where  an  incoming  flow  of  Mach  5  and  a  unit  Reynolds 
number  of  37xl0^/m  is  colliding  with  an  inclined  wall  (shock  generator),  generating  a  shock  that  impinges  on  the 
bottom  flat  plate  interacting  with  the  flat  plate  boundary  layer.  The  experimental  set-up  and  the  basic  flow  features 
in  the  vicinity  of  the  SWBEI  are  shown  in  Figure  1. 


Figure  1:  Experimental  set-up  of  Schiilein  (left)  (adapted  from  Schiilein^)  and  schematics  of  the  flow  physics 

near  the  SWTBLI  region  (right) 

The  experimental  data  include  wall  pressure,  skin  friction  coefficient  and  velocity  profiles  for  three  different 
shock  generator  angles:  6,  10  and  14  degrees.  The  cases  for  angles  10  and  14  degrees  correspond  to  higher  shock 
intensity  levels  compared  with  case  6,  where  the  boundary  layer  is  fully  separated. 
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Figure  2:  Fine  grid  at  14  degrees  shock-generator  angle  for  the  SWTBLI. 


Figure  3:  Mach  number  contours  for  the  14  degrees  shock-generator  angles  results  shown  for  the  WENO  3*^^ 

order  scheme. 

Computations  were  performed  for  shock-generator  angles  10  and  14  degrees  using  the  MUSCL  and  WENO 
schemes.  The  time  integration  was  performed  by  a  third-order  explicit  Runge-Kutta  scheme.  In  the  framework  of 
RANS,  a  grid  convergence  study  was  performed  using  coarse,  medium  and  fine  grids  consisting  of  59,405,  96,879 
and  179,997  elements,  respectively.  The  grids  are  locally  refined  to  capture  the  boundary  layer  effects  and  the 
impinging  shock;  Figure  2  illustrates  the  fine  grid  for  case  14;  the  computed  Mach  number  contour  levels  are  shown 
in  Figure  3  for  the  3^^- order  WENO  scheme. 

The  wall  pressure  and  skin  friction  coefficient  on  the  lower  wall  (flat  plate)  for  the  2"^-order  MUSCF  (M2)  and 
3'^^-order  WENO  (W3)  schemes  are  compared  with  the  experimental  data  for  the  14  degrees  case  on  the  fine  grid 
(Figure  4).  Overall  the  agreement  is  acceptable  for  the  pressure  loads  apart  from  the  under-prediction  near  the 
separation  point.  In  terms  of  skin-friction,  considerable  under-predictions  are  present  downstream  of  the  separation 
to  the  end  of  the  flat  plate.  The  W3  solution  recovers  better  featuring  smaller  oscillations  in  the  recirculation  region 
(300<x<350). 

Two  different  techniques  were  employed  for  the  discretization  of  the  turbulence  transport  equation:  i)  using  the 
same  reconstruction  as  for  the  mean  flow  equations'^  (labeled  as  coupled  approach);  ii)  using  a  first  order  upwind 
method  (labeled  as  decoupled  approach).  The  results  are  shown  in  Figure  5,  where  the  coefficient  of  skin  friction  is 
plotted  for  the  10-degree  case  on  the  medium  grid.^^  Furthermore,  the  viscous  gradient  computation  method  has 
been  investigated  in  conjunction  with  the  Green  Gauss  (GG)  and  the  least-square  method  (ESQ).  It  was  found  that 
the  profile  obtained  by  the  ESQ  method  provides  better  agreement  with  the  experiment  near  the  recirculation  region. 
The  Mach  number  boundary  layer  profile  is  plotted  for  the  10  degrees  case  at  two  stations  in  Figure  6  showing 
satisfactory  agreement  with  the  experiment. 
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Figure  4:  Dimensionless  wall  pressure  on  the  lower  wall  (left)  and  skin  friction  coefficient  (right).  The 
computations  using  the  MUSCL  l^^^-order  (M2)  and  WENO  3*^‘^-order  (W3)  schemes  are  compared  with  the 

experiment. 


Figure  5:  Skin  friction  coefficient  on  the  lower  wall  for  the  10  degrees  case,  for  the  M3  and  W3  schemes  with 
two  different  viscous  gradient  methods,  the  Green  Gauss  (GG)  and  the  Least  Square  method  (LSQ);  coupled 
results  are  shown  on  the  left  and  decoupled  results  are  shown  on  the  right. 
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Figure  6:  Mach  number  profiles  for  the  10  degrees  case  at  two  different  stations;  left  at  x=396mm  and  right 

at  x=426mm. 


B.  Compression  ramp 

ILES  has  been  performed  to  simulate  the  shock  wave/turbulent  boundary  layer  interaction  (SWTBLI)  over  a 
compression  ramp.  The  structured  finite  volume  code  CNS3D  will  be  used,  which  numerical  background  is  given  in 
section  §III  B.  The  case  is  based  on  the  DNS  study  of  Wu  and  MarW^  where  an  incoming  flow  of  Mach  2.9  and  a 
Reynolds  number  of  38.7x10^  based  on  the  freestream  properties  and  boundary  layer  height  (^5)  collides  with  a 
surface  as  the  wall  in  inclined  abruptly  by  24°,  generating  a  shock  that  interacts  with  the  incoming  turbulent  flow. 
The  boundary  layer  forms  a  separation  bubble  with  its  size  being  dictated  by  the  intensity  of  the  incoming  turbulent 
flow  and  the  strength  of  the  formed  shockwave. 

The  effects  of  grid  resolution  was  investigated  using  three  different  grids  (see  Table  1).  Each  level  of  grid 
refinement  used  a  smaller  y+  value  of  the  first  point  from  the  wall.  Previous  investigations^^  regarding  low-Mach 
compressible  turbulent  boundary  layers,  have  shown  that  y+  values  of  around  2  is  sufficient,  at  least  for  the  9^^-order 
WENO  scheme.  Furthermore,  the  5^^-order  MUSCE  in  conjunction  with  the  low-Mach  number  correction  of^  is 
expected  to  provide  accurate  results.  The  low-Mach  correction  is  required  for  capturing  the  subsonic  region  of  the 
incoming  turbulent  boundary  layer  and  its  influence  decreases  linearly  to  zero  as  the  transonic  boundary  limit  is 
reached. 


Nx 

Nz 

Ny 

z+ 

Coarse 

600 

96 

72 

2 

Medium 

840 

120 

96 

1 

Fine 

1128 

168 

120 

0.5 

DNS^^ 

1024 

160 

128 

0.2 

Table  1:  Grid  resolution  in  number  of  cells 


The  size  of  the  computational  domain  used  is  compared  to  that  of  the  DNS  in  Table  2.  Note  that  the  height  of  the 
domain  at  the  inflow  location  has  remained  the  same.  However  a  larger  streamwise  length  was  utilized  in  order  to 
accommodate  the  spatially  developing  turbulent  boundary  layer. 
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L 

L 

L 

X 

y 

z 

ILES 

21.4 

3 

5 

DNS^^ 

15.4 

2.2 

5 

Table  2:  Size  of  computational  domain  (in  terms  of  <5). 


Periodic  boundary  conditions  were  used  in  the  spanwise  (y)  direction.  In  the  wall-normal  (z)  direction,  a  no-slip 
isothermal  wall  (Tw=309K)  was  used.  Supersonic  outflow  was  prescribed  at  the  outlet,  while  free  stream  conditions 
(far-field)  were  assigned  at  the  upper  boundary  opposite  to  the  wall.  The  boundary  condition  at  the  inlet  requires  of 
accurately  assigning  a  turbulent  boundary  layer.  A  synthetic  turbulent  digital  filter  approach^^  was  further  developed 
in  conjunction  with  the  present  numerical  framework^^  to  generate  incoming  turbulent  boundary  layer  data.  It  is 
shown  to  work  satisfactorily  and  has  since  been  used  successfully  in  a  number  of  different  test-cases^^’^^’"^^.  Note  that 
since  the  numerical  scheme  still  needs  to  capture  and  resolve  the  synthetic  inflow  perturbations  that  lead  to  a 
turbulent  flow,  the  length  of  the  upstream  domain  is  increased  by  approximately  five  boundary  layer  heights. 


Figure  7:  Computational  grids  for  the  ramp  case. 
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The  computational  grids  are  shown  in  Figure  7.  The  coarse  grid  comprises  4,147,200  cells,  while  the  medium 
and  fine  grids  comprise  of  9,676,800  and  22,740,480,  respectively. 

In  order  to  reduce  the  computational  cost,  particularly  on  the  fine  grid,  a  variable  limiter  approach  was 
employed.  The  9^^-order  WENO  is  more  expensive  than  the  2”^  and  5^^-order  schemes  due  to  the  increased  number 
of  candidate  stencils.  Using  a  reduced  order  of  accuracy  in  the  outer  region  of  the  supersonic  laminar  flow  does  not 
negatively  impact  on  the  accuracy  of  the  results  in  the  near  wall  region.  Any  disturbances  that  arise  from  the  inner 
turbulent  boundary  layer  are  steadily  advected,  but  more  rapidly  dissipated.  A  2”^-order  MUSCL  limiter  is  used 
(Monotonized  Central,  MC)  in  the  outer  wall  regions  of  the  flow  allowing  for  the  more  computationally  expensive 
5^^-order  MUSCL  and  9^^-order  WENO  schemes  to  be  used  in  the  near  wall  region.  The  computational  cells  located 
further  than  two  boundary  layer  thicknesses  (2S)  from  the  closest  wall  point  are  deemed  marginal  and  assigned  to 
the  lower  numerical  accuracy  region.  Error!  Reference  source  not  found,  shows  the  resulting  regions  of  low  and 
high-order  of  accuracy  as  applied  to  the  compression  ramp  geometry  and  settings  described. 


Red  Region:  MUSCL  5*^  or  WENO  9*^  order  limiters 
z/6<2 


Figure  8:  Variable  limiter  regions 

The  reduction  of  computational  cost  gained  by  implementing  variable  limiter  regions  depends  on  the  block 
partitioning  of  the  domain  conducted  for  parallel  processing.  In  the  case  of  the  9^^-order  WENO  the  variable  limiter 
region  approach  was  found  to  speed  up  the  computations  by  approximately  10-20%. 

As  remarked  earlier,  the  size  of  the  separation  bubble  will  highly  depend  on  the  nature  of  the  turbulent  flow 
resolved.  This  is  because  the  turbulent  boundary  layer  will  frequently  deposit  large  amounts  of  momentum  near  the 
wall  surface,  which  acts  to  push  the  “growing”  separation  bubble  back.  In  the  presence  of  a  laminar  boundary  layer, 
a  separation  bubble  is  expected  to  gradually  increase  in  size  by  gradually  “creeping”  upstream  via  the  subsonic  part 
of  the  boundary  layer.  Turbulence  acts  to  impede  on  this  process  by  moving  supersonic  flow  of  high  momentum  and 
kinetic  energy  closer  to  the  wall  surface. 

Figure  9  gives  a  qualitative  comparison  of  the  level  of  turbulence  resolved  by  each  numerical  method  and  grid 
resolution  examined  by  plotting  the  iso-surfaces  of  the  Q-criterion.  As  expected,  each  level  of  grid  refinement  leads 
to  an  increase  in  the  amount  and  strength  of  vortex  cores  resolved.  Noticeably,  the  9^^-order  WENO  on  the  coarse 
grid  appears  capable  of  resolving  a  similar  number  of  vortices  in  the  region  preceding  the  separation  bubble  as  the 
5^^-order  MUSCL.  An  integral  length  scale  of  2Ax  in  the  streamwise  direction  was  chosen  for  the  digital  filter 
turbulent  inflow  technique.  Past  the  shockwave  front,  the  grid  refinement  appears  to  have  little  to  no  effect  on  the 
structure  of  the  vortices  resolved  by  the  9^^-order  WENO.  On  the  contrary,  the  5^^-order  MUSCL  shows  a  gradual 
improvement  in  the  resolved  turbulent  structure,  which  eventually  resembles  that  obtained  by  the  9^^-order  WENO. 
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MUSCL  5™  +  LM  WENO  9™ 


b)  Medium 


c)  Fine 


Figure  9:  Iso-surfaces  of  ^-criterion;  2*  (<5/Uoo)^=2  &  5  colored  by  />/>9ooe[-0.2,2]  as  xz-plane  contour  plot. 
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MUSCL  5™  +  LM 


a)  Coarse 


b)  Medium 


c)  Fine 

Figure  10:  Iso-surfaces  depicting  shock-bubble  interaction;  blue:  separation  bubble  />u=0,  orange:  shock 
pressure  P/Pw~1.765  &  2.353,  xz-plane:  contrours  of/>/>9ooe[0.4,3]. 


WENO  9™ 
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Figure  10  illustrates  the  interaction  of  the  separation  bubble  and  shockwave.  The  most  important  factor  that 
determines  the  physical  processes  occurring  here  is  the  incoming  turbulent  flow,  which  relies  largely  on  the 
numerical  scheme  used  as  already  discussed.  The  “streaky”  appearance  of  the  separation  bubble’s  leading  edge  is  a 
result  of  the  incoming  turbulent  flow,  which  deposits  energy  from  the  freestream  onto  the  near-wall  surface.  The 
effect  is  to  push  random  areas  of  the  bubble  downstream,  thus  causing  the  observed  streaks.  The  shockwave  pressure 
front  is  found  to  be  more  unsteady-like  (as  evident  by  the  wavy-like  structure)  in  the  frontal  region  of  the  separation 
bubble.  However  this  becomes  less  evident  as  the  grid  is  refined.  The  cause  for  this  unexpected  result  is  the 
turbulent  integral  length  scales  produced  by  the  turbulent  inflow  digital  filter  technique,  which  relies  on  the  inflow 
cell  length  to  give  a  value  to  the  streamwise  turbulent  integral  length  scale.  Increasing  the  number  of  computational 
cells  leads  to  estimating  smaller  turbulent  integral  length  scales  fed  into  the  digital  filter.  This  leads  inevitably  to 
smaller  “large”  turbulent  integral  scales  produced  in  the  flow,  which  cannot  “carry”  as  much  energy  from  the 
freestream  to  the  near-wall  region,  thus  greatly  affecting  the  strong  shockwave  pressure  front  formed. 


Figure  11:  Contours  of  time-averaged  TKE  puu  !  2pJJl^ ;  (left:  ILES  coarse  grid  W9,  right:  DNS^^). 


The  time-averaged  turbulent  kinetic  energy  (TKE)  resolved  by  the  9^^-order  WENO  scheme  on  the  coarse  grid, 
as  seen  in  Figure  11,  is  consistent  to  that  obtained  by  the  DNS  of  Wu  and  Martin^^,  albeit  the  maximum  TKE  is 
found  over  a  lesser  area.  Nevertheless,  this  is  still  a  very  encouraging  result  considering  the  difference  in  the  total 
number  of  cells  used  in  either  case. 


e 


6 
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Figure  12:  Location  of  velocity  profiles  used  for  comparisons  with  DNS. 
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Figure  12  shows  the  locations  along  which  time-averaged  velocity  profiles  are  compared  to  DNS.  Three 
locations  are  chosen:  i)  the  first  is  located  in  the  incoming  turbulent  inflow  and  acts  as  a  good  measure  of  the  quality 
of  the  resolved  mean  turbulent  profile;  ii)  the  next  is  located  prior  to  the  compression  corner  but  within  the  leading 
front  of  the  separation  bubble  and  acts  as  an  indicator  of  how  well  the  separation  bubble  is  captured;  and,  finally,  iii) 
the  last  is  located  just  prior  to  the  domain  exit  (outflow  boundary  condition)  and  reveals  whether  the  correct 
turbulent  profile  and  shock  position  are  captured. 

The  mean  turbulent  velocity  profiles  (Figure  13)  show  a  satisfactory  agreement  to  DNS.  Generally,  the  9^^-order 
WENO  captures  a  sharper  velocity  gradient  in  the  near-wall  region  than  the  5^^-order  MUSCL  scheme.  Most 
interestingly,  however,  the  9^^-order  WENO  results  obtained  on  the  coarse  grid  agree  best  to  the  DNS.  This  is  due  to 
the  dependence  of  the  turbulent  integral  length  scale  used  by  the  digital  filter  on  the  grid  resolution,  as  previously 
explained. 
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The  impact  on  the  results  of  the  decreasing  turbulent  length  scales  associated  with  the  increased  grid  resolution 
becomes  more  evident  in  Figure  14.  The  5  -order  MUSCL  seems  to  be  less  sensitive  to  the  grid  resolution, 
however,  this  is  partly  due  to  its  reduced  ability  to  resolve  features  of  two  cell  lengths  imposed  at  the  inflow  by  the 
digital  filter.  In  the  post-shock  region  and  just  prior  to  the  outflow.  Figure  15  shows  that  the  5^^-order  MUSCL 
scheme  has  given  a  better  agreement  to  the  DNS  than  the  9^^-order  WENO.  It  is  difficult  to  ascertain  a  reason  for 
this  occurrence  at  this  early  stage,  but  the  erroneous  results  obtained  upstream  could  be  partly  to  blame.  On  another 
note,  the  transition  from  the  high-order  to  the  low-order  region  does  not  appear  to  produce  any  kind  of  unphysical 
behavior  at  least  to  the  mean  velocity  profile  examined  in  Figure  15. 
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Pw^*.  P«/P. 


The  mean  wall  pressure  distribution  obtained  by  the  5^^-order  MUSCL  with  low-Mach  correction  and  the  9^^- 
order  WENO  are  compared  to  DNS^^  and  experiment"^^  in  Figures  17  and  18,  respectively.  Both  schemes  give  the 
best  agreement  on  the  coarse  grid,  as  evident  by  the  pressure  profile  at  the  separation  bubble  that  begins  at  x/S=-5  up 
to  the  corner  at  x/S  =0.  In  the  remaining  regions,  both  schemes  on  all  grids  give  wall  pressure  values  that  fall  within 
those  obtained  by  DNS  and  experiment. 
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V.  Conclusions 


The  capability  of  high-order,  high-resolution  methods  is  examined  in  the  context  of  shock  wave/turbulent 
boundary  layer  interaction  (SWTBLI)  from  moderate  to  high  Reynolds  numbers.  An  unstructured  solver  was  used  to 
model  an  impinging  shock  at  a  sufficiently  high  Reynolds  numbers  that  required  application  of  the  Spalart-Allmaras 
(SA)^^  turbulence  model.  Overall  results  are  found  to  be  encouraging,  albeit  a  number  of  issues  regarding  numerical 
stability  of  the  simulations  required  investigation.  Agreement  of  the  pressure  loads  is  found  to  be  acceptable  apart 
from  the  under-prediction  near  the  separation  point.  In  terms  of  skin-friction  considerable  under-predictions  are 
present  downstream  of  the  separation  to  the  end  of  the  flat  plate.  The  B'^^-order  WENO  solution  recovered  better 
featuring  smaller  oscillations  in  the  recirculation  region  (300<x<350).  It  is  also  found  that  the  profile  obtained  by  the 
least-square  method  (ESQ)  provides  better  agreement  with  the  experiment  near  the  recirculation  region.  The  Mach 
number  boundary  layer  profile  also  showed  satisfactory  agreement  with  the  experiment.  With  regards  to  the 
turbulence  model,  coupling  it  to  the  same  numerical  schemes  used  during  the  estimation  of  the  inviscid  fluxes, 
namely  the  reconstruction  method,  resulted  in  superior  accuracy  compare  to  the  traditionally  used  2”^-order  central 
method. 

The  accuracy  of  the  5^^-order  MUSCE,  in  conjunction  with  the  low-Mach  correction,  and  the  9^^-order  WENO 
schemes  were  investigated  in  SWTBEI  over  a  compression  ramp.  The  Reynolds  number  was  sufficiently  low  to 
allow  for  lEES  to  be  conducted.  The  results  showed  a  great  deal  of  sensitivity  to  the  resolved  incoming  turbulent 
flow  produced  using  the  digital  filter  technique.  Time-averaged  streamwise  velocity  and  wall  pressure  profiles 
suggest  that  the  finer  (smaller)  largest  integral  scales  produced  by  the  digital  filter,  as  a  result  of  the  dependence  of 
the  estimation  of  the  turbulent  length  scales  to  the  inflow  cell  size,  produces  a  “weaker”  structure  that  can  move  less 
kinetic  energy  from  the  freestream  to  the  inner-wall  regions  of  the  boundary  layer.  As  a  result,  increasing  the  grid 
resolution  causes  the  leading  edge  of  the  separation  bubble  to  move  un-physically  further  upstream.  The  5^^-order 
MUSCE  was  not  accurate  enough  to  sufficiently  resolve  the  perturbations  produced  by  the  digital  filter  based  on 
mean  turbulent  integral  scales  of  twice  the  cell  length.  Finally,  a  variable  limiter  region  approach  was  implemented 
that  reduced  the  computational  requirements  by  10%  to  20%  depending  on  the  numerical  scheme  and  a  careful 
(MPI)  block  decomposition. 
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